Computationally Efficient Finite Impulse Response Comb Filtering

ABSTRACT

A method and system that remove an unwanted signal and its harmonics from an input signal in a computationally efficient manner are disclosed. Embodiments include processing the FFT matrix to selectively zero-out rows of the matrix before multiplying the matrix with the Inverse FFT (IFFT) matrix. The resulting product (which is a sparse matrix) is then used to generate coefficients for a linear Finite Impulse Response (FIR) filter to process the input. The filtered output signal has the unwanted signal and its harmonics removed with minimal effect on a desired signal. The method produces a stable, physically realizable filter, requiring fewer computations than current methods.

STATEMENT REGARDING FEDERALLY-SPONSORED RESEARCH AND DEVELOPMENT

Statement under M.P.E.P. §310. The U.S. government has a paid-up licensein this invention and the right in limited circumstances to require thepatent owner to license others on reasonable terms as provided for bythe terms of contract W15P7T-12-C-F600 awarded by the United States AirForce.

FIELD OF THE INVENTION

Embodiments described herein generally relate to digital filtering.

BACKGROUND OF THE INVENTION

Comb filters have use in various applications. One application includesthe removal of a fundamental sinusoidal signal and its harmonics from asignal of interest (SOI).

Existing comb filtering techniques have practical limitations. A firsttechnique includes the use of a finite impulse response (FIR) filter.Although this technique is fast and simple, it typically has a poorfrequency response. A second technique involves the use of an infiniteimpulse response (IIR) filter. Theoretically, a sharp notch filter canbe created with a lower-order IIR filter than produced with an FIRfilter. However, the poles of such a filter must be very close to theunit circle, which causes stability and performance problems when thefilter coefficients are quantized for implementation. A third techniqueinvolves transforming a block of data samples into the frequency domainusing a Fast Fourier Transform (FFT) (or, more generally, a discreteFourier Transform), and zeroing the bins (“zero-binning”) of thefrequencies associated with the interference and its harmonics. Aninverse transformation back to the time domain produces the filteredoutput data. This technique produces acceptable results, but has limitedapplication due of its computational complexity.

SUMMARY OF THE INVENTION

Embodiments for FFT comb filtering are provided. Embodiments operate byzeroing every M^(th) row of an N-point FFT matrix and multiplying theresult by an IFFT matrix before operation on the input data. If N/M isan integer, multiplying the modified FFT matrix with the IFFT matrixproduces a sparse matrix that can be manipulated to form a FIR filter.The input data is passed sequentially through the FIR filter, and thespectral contents of a set of uniformly spaced frequencies are nulledand removed from the data. This novel approach accomplishes thefiltering without having to transform the input data into the frequencydomain.

BRIEF DESCRIPTION OF THE DRAWINGS/FIGURES

The accompanying drawings, which are incorporated herein and form a partof the specification, illustrate the present invention and, togetherwith the description, further serve to explain the principles of theinvention and to enable a person skilled in the pertinent art to makeand use the invention.

FIG. 1 is a block diagram that illustrates a FFT-based method of combfiltering.

FIG. 2 is a block diagram that illustrates an example comb filteraccording to an embodiment.

FIG. 3 is a block diagram of another example comb filter according to anembodiment.

FIG. 4 is a flow chart that illustrates a comb filtering processaccording to an embodiment.

DETAILED DESCRIPTION OF THE INVENTION

FIG. 1 is a block diagram that illustrates a FFT-based comb filter 100.As shown in FIG. 1, comb filter 100 includes a Fast Fourier Transform(FFT) module 102, a zero-binning module 104, and an Inverse FFT (IFFT)module 106. Comb filter 100 may be used to remove an interferingcontinuous wave (CW) signal.

For purposes of this discussion, the term “module” shall be understoodto include at least one of software, firmware, and hardware (such as oneor more circuits, microchips, or devices, or any combination thereof),and any combination thereof. In addition, it will be understood thateach module can include one, or more than one, component within anactual device, and each component that forms a part of the describedmodule can function either cooperatively or independently of any othercomponent forming a part of the module. Conversely, multiple modulesdescribed herein can represent a single component within an actualdevice.

Comb filtering by filter 100 begins with FFT module 102 receiving a timedomain input signal 108. In an embodiment, input signal 108 includes adesired signal and an interfering CW signal. In an embodiment, inputsignal 108 includes an N-point block of input samples. FFT module 102acts on input signal 108 to produce a frequency domain output signal110.

In an embodiment, output signal 110 represents an N-point FFT of inputsignal 108. Generally, an N-point FFT module with a sample frequency,f_(s), will have N frequency bins that contain frequency componentsspaced at intervals of:

$\begin{matrix}{f_{d} = {\frac{f_{s}}{N}.}} & (1)\end{matrix}$

When the CW frequency and its harmonics fall exactly at integer, k,multiples of f_(d), these frequency components, kf_(d), will becompletely contained in those bins.

Output signal 110 is fed into zero binning module 104. Zero binningmodule 104 acts on the bins of output signal 110 that include theinterfering CW frequency and its harmonics. Specifically, zero binningmodule 104 zeroes the bins that include the CW frequency and itsharmonics to produce signal 112. When these bins are zeroed, theinterfering CW signal is eliminated.

Signal 112 is then provided to IFFT module 106. IFFT module 106transforms signal 112 from the frequency domain to the time domain byapplying an N-point IFFT to generate a filtered output signal 114.

While the conventional comb filtering technique described in FIG. 1produces acceptable results generally, it suffers from excessivecomputational complexity. Specifically, the technique requirescollecting N samples of data, performing an FFT on this data, zeroingevery M^(th) value of the resulting FFT, and then performing an inverseFFT (IFFT) to obtain the time domain filtered result. The FFT and IFFTeach requires on the order of N log₂N multiplications, resulting in 2Nlog₂N multiplications for the entire process. As a result, practicalhardware limitations (due to the size and speed of FFT manipulationsthat are required) eliminate this technique from usage for all butlow-bandwidth applications.

Embodiments enable a comb filter that does not suffer from deficienciesof existing comb filtering techniques. Specifically, embodimentsrecognize that FFT comb filtering can be performed by applying an inputsignal through a module that implements the multiplication of the inputsignal with the product of an inverse FFT matrix and appropriatelyzeroed FFT matrix. Specifically, for an N-sample input signal, FFT combfiltering can be realized by multiplying the input signal with theproduct of an inverse FFT matrix and an FFT matrix, with every M^(th)row of the FFT matrix set to zero, with the condition that N and M areselected such that N/M is an integer. As such, comb filtering of theinput signal is realized without performing a FFT or an IFFT on theinput signal.

A mathematical explanation of an embodiment is provided below for thepurpose of illustration. As would be understood by a person of skill inthe art, embodiments are not limited by this mathematical explanation.

Let A and A⁻¹ be the size N FFT and IFFT matrices. A and A⁻¹ may be thematrices implemented by FFT module 102 and IFFT module 106, for example.A is given by

$\begin{matrix}{{\frac{1}{\sqrt{N}}\begin{bmatrix}1 & 1 & \ldots & 1 \\1 & W_{N}^{1} & \ldots & W_{N}^{({N - 1})} \\1 & W_{N}^{2} & \ldots & W_{N}^{2{({N - 1})}} \\\vdots & \vdots & \vdots & \vdots \\1 & W_{N}^{({N - 1})} & \ldots & W_{N}^{{({N - 1})}{({N - 1})}}\end{bmatrix}},{{{where}\mspace{14mu} W_{N}} = {^{{{- j}\; \frac{2\pi}{N}}\;}.}}} & (2)\end{matrix}$

Embodiments recognize that zeroing every M^(th) value of the product ofthe FFT matrix with the input vector, x, is equivalent to multiplyingthe input vector with the FFT matrix with every M^(th) row of the FFTmatrix zeroed. This is further described below.

Let A_(r) be the FFT matrix with every M^(th) row zeroed. Thus, A_(r) isgiven by:

$\begin{matrix}{{\frac{1}{\sqrt{N}}\begin{bmatrix}0 & 0 & \ldots & \; & 0 \\1 & W_{N}^{1} & \ldots & \; & W_{N}^{({N - 1})} \\1 & W_{N}^{2} & \; & \; & W_{N}^{2{({N - 1})}} \\\vdots & \vdots & \; & \; & \vdots \\0 & 0 & \ldots & \; & 0 \\\vdots & \vdots & \vdots & \; & \vdots \\\; & \; & \; & \; & \; \\1 & W_{N}^{({N - 1})} & \ldots & \; & W_{N}^{{({N - 1})}{({N - 1})}}\end{bmatrix}}.} & (3)\end{matrix}$

Let C_(r) be the result of the multiplication of A⁻¹ and A_(r):

C _(r) =A ⁻¹ A _(r)  (4)

A_(r) can also be represented as the matrix (A−A_(e)), where A_(e) isthe matrix A_(e) with all its rows zeroed out except the first row andevery M^(th) row. Then (4) can also be expressed as:

$\begin{matrix}\begin{matrix}{C_{r} = {A^{- 1}\left( {A - A_{e}} \right)}} \\{= {I - {A^{- 1}{A_{e}.}}}}\end{matrix} & (5)\end{matrix}$

Since only every M^(th) element of the columns of A_(e) is non-zero, thesummation is only over 0 to (N/M−1). The elements of matrix if A⁻¹A_(e)are sums of N/M products and can be expressed as:

$\begin{matrix}{{A^{- 1}{A_{e}({ij})}} = {\frac{1}{N}{\sum\limits_{n = 0}^{\frac{N}{M} - 1}{W_{N}^{inM}{W_{N}^{- {jnM}}.}}}}} & (6)\end{matrix}$

This simplifies to:

$\begin{matrix}{{A^{- 1}{A_{e}({ij})}} = {\frac{1}{N}{\sum\limits_{n = 0}^{\frac{N}{M} - 1}{W_{N/M}^{{({i - j})}n}.}}}} & (7)\end{matrix}$

Thus, A⁻¹A_(e) can be written as:

$\begin{matrix}{{A^{- 1}A_{e}} = \left\{ \begin{matrix}\frac{1}{M} & {{{for}\mspace{14mu} n} = {{{mod}\left( {{i - j},\frac{N}{M}} \right)} = 0}} \\0 & {{otherwise}.}\end{matrix} \right.} & (8)\end{matrix}$

Then, the ij^(th) element of C_(r) can be written as:

$\begin{matrix}{{C_{r}({ij})} = \left\{ \begin{matrix}{1 - \frac{1}{M}} & {i = j} \\{- \frac{1}{M}} & {{n = {{{mod}\left( {{i - j},\frac{N}{M}} \right)} = 0}},{{{and}\mspace{14mu} i} \neq j}} \\0 & {{otherwise}.}\end{matrix} \right.} & (9)\end{matrix}$

Thus, C_(r) is a multi-diagonal matrix with 2M−1 diagonals, includingthe main diagonal. The diagonals are spaced every N/M columns or rowsapart. The main diagonal elements have the value (1−1/M) and theoff-diagonal elements with non-zero elements all have the value −1/M.

Accordingly, if the input vector, x, is multiplied by the matrix C_(r),the first value of the resulting vector can be written as:

$\begin{matrix}{{\hat{x}}_{0} = {x_{0} - {\frac{1}{M}{\sum\limits_{k = 0}^{M - 1}{x_{{- k}\; \frac{N}{M}}.}}}}} & (10)\end{matrix}$

In general, the k^(th) value of the resulting vector can be written as:

$\begin{matrix}{{{\hat{x}}_{k} = {x_{k} - {\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}x_{k - {j\; \frac{N}{M}}}}}}},} & (11)\end{matrix}$

for k equal 0 to N/M−1.

Because of the cyclic nature of the matrix, C_(r), only the first N/Mvalues need to be calculated as in (11). The rest of the (N−N/M) valuescan be calculated using the first N/M values as:

$\begin{matrix}{{{\hat{x}}_{k} = {{\hat{x}}_{{mod}{({k,{N/M}})}} + x_{k} - x_{{mod}{({k,{N/M}})}}}},} & (12)\end{matrix}$

for k equal N/M to N−1.

As shown above, the total number of multiplications for the abovedescribed embodiment reduces to N/M, and all other operations areadditions. Thus, embodiments provide an equivalent method ofimplementing comb filtering which is a factor of 2M log₂/N moreefficient than the conventional FFT method.

As a result, computation time can be reduced dramatically usingembodiments. For example, for low bandwidth systems, N is typicallyequal to 4096 and M is typically equal to 32. With these values, thefactor of improvement in computational complexity equals 768. For largerbandwidth systems, N may be 2¹⁹ or 524288. Embodiments reduce thecomputations by a factor of 1216 less multiplications. Further,embodiments eliminate the requirement that a radix 2 value be used forthe size of the FFT. This allows greater flexibility in implementation.

There is still the requirement that all N samples must be collectedbefore the outputs can be obtained. However, examination of (11) showsthat the filtering operation is reduced to a linear FIR filter. This ishighlighted in (13) as

$\begin{matrix}{{{\hat{x}}_{0} = {x_{0} - {\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}x_{k - {j\; \frac{N}{M}}}}}}}{y_{k} = {x_{k} - {\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}{x_{k - {j\; \frac{N}{M}}}.}}}}}} & (13)\end{matrix}$

The comb filtered output can be obtained by filtering the input with thecoefficients of the FIR filter in (13). While the filter operates over awide range of N input values, only M values are involved in each outputcalculation. Only one multiplication and M additions are required. If Mis radix 2, no multiplications, only bit shifts, are required. Once Ninput samples have been obtained, i.e., filter startup has completed,there is no delay in getting an output value with each new input sample.

FIG. 2 is a block diagram that illustrates an example comb filter 200according to an embodiment. Example comb filter 200 is provided for thepurpose of illustration only and is not limiting of embodiments. Asshown in FIG. 2, example comb filter 200 includes a N-stage shiftregister 202, an adder module 204, a multiplier module 206, and an addermodule 208.

N-stage shift register 202 is configured to receive an input signal 210.In an embodiment, input signal 210 includes an N-point block of inputsamples, which are input into N-stage shift register 202 in a serialmanner. N-shift register 202 is configured to produce an output signal212, which includes data values from select registers of N-shiftregister 202. In an embodiment, output signal 212 includes M valuesselected such that N/M is an integer.

Output signal 212 is provided to adder module 204, which is configuredto add the data values contained in output signal 212 to generate asignal 216. Signal 216 is then provided to multiplier module 206, whichis configured to multiply signal 216 by a scalar to produce a signal218.

Subsequently, signal 218 is provided to adder module 208, which addssignal 218 to a current sample (time index k) of input signal 210 togenerate output signal 220. Output signal 220 represents the outputsample at time index k.

FIG. 3 is a block diagram of another example comb filter 300 accordingto an embodiment. As shown in FIG. 3, example comb filter 300 includes aplurality of registers 302, an adder module 304, a multiplier module306, and an adder module 308.

Registers 302 are configured to receive and store sequential samples ofan input signal 310. In an embodiment, registers 302 includes[(M−1)N/M+1] registers, with N/M being an integer.

Adder module 304 is configured to sum values from select registers ofthe plurality of registers 302 to generate a signal 312. Signal 312 ismultiplied by a scalar via multiplier module 306 to produce a signal314. In an embodiment, the scalar is equal to −1/M, where M is aninteger selected such that N/M is an integer.

Adder module 308 is configured to sum signal 314 and a current sample(time index k) of input signal 310 to form a current sample (time indexk) of a filtered output signal 316 of comb filter 300.

FIG. 4 is a flow chart that illustrates a comb filtering process 400according to an embodiment. Process 400 is provided for the purpose ofillustration and is not limiting. As shown in FIG. 4, process 400 beginsat the start (402), and then proceeds to the designer's selection ofthree parameters, M, f_(s), and N (404), which determine thecharacteristics of the filter. The factor M roughly equates to the Qfactor of a conventional filter design. The choice of M (for a specificsampling frequency and number of samples) determines the width of theindividual filter notches. N/M determines the number of notches.Additionally, the sampling frequency, f_(s), and the number of points,N, that would be used for an equivalent FFT-IFFT process are selected.The shift register is then configured according to the selection of theprevious parameters and the filtering equation, (13), determined by thepresent invention (406). Then, the input data values are input to thefilter (408), and after processing are available at the output (410).The comb filtering process is completed at the end (412).

What is claimed is:
 1. A method for filtering an input signal,comprising: receiving an input signal vector that corresponds to theinput signal; and applying the input signal vector to a filter thatemulates a product of an inverse Fast Fourier Transform (IFFT) matrixand a selectively zeroed Fast Fourier Transform (FFT) matrix, whereinthe selectively zeroed FFT matrix has select rows set to zero to reducean interfering signal contained in the input signal.
 2. The method ofclaim 1, wherein the input signal vector includes N samples of the inputsignal, and wherein every M^(th) row of the selectively zeroed FFTmatrix is set to zero, and wherein N/M is an integer.
 3. The method ofclaim 1, wherein the select bins include bins that include components ofa fundamental frequency of the interfering signal.
 4. The method ofclaim 1, wherein the select bins includes bins that include harmoniccomponents of the interfering signal.
 5. The method of claim 1, whereinapplying the input signal vector to the filter comprises applying theinput signal to a filter that implements the following equation:$y_{k} = {x_{k} - {\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}x_{k - {j\; \frac{N}{M}}}}}}$where y_(k) represents a filtered output sample at index k, x_(k)represents a sample of the input signal at index k, and M is an integer.6. The method of claim 5, wherein M is selected to notch a fundamentalfrequency and harmonic components of the interfering signal.
 7. A combfilter, comprising: a shift register comprising a plurality ofregisters, each configured to receive a respective sample of an inputsignal; a first adder module configured to sum input sample values fromselect registers of the plurality of registers to generate a firstsignal; a multiplier module configured to multiply the first signal by ascalar to generate a second signal; and a second adder module configuredto sum the second signal with a current sample value of the input signalto generate a current sample value of a filtered output signal.
 8. Thecomb filter of claim 7, wherein the input signal includes N samples, andwherein the plurality of registers include (M−1)(N/M)+1 registers, andwherein N and M are integers and N/M is an integer.
 9. The comb filterof claim 7, wherein the scalar is equal to −1/M, wherein M is such thatN/M is an integer, and wherein N represents a number of samples of theinput signal.
 10. The comb filter of claim 7, wherein the scalar isselected based on a product of an inverse Fast Fourier Transform (IFFT)matrix and a selectively zeroed Fast Fourier Transform (FFT) matrix. 11.The comb filter of claim 7, wherein the select registers are chosenbased on a product of an inverse Fast Fourier Transform (IFFT) matrixand a selectively zeroed Fast Fourier Transform (FFT) matrix.
 12. Thecomb filter of claim 7, wherein the shift register is implemented usinga Field Programmable Gate Array (FPGA).
 13. A method of comb filteringan input signal, comprising: receiving an input signal vector thatcorresponds to the input signal; and applying the input signal vector toa linear Finite Impulse Response (FIR) filter, wherein an output signalvector of the FIR filter is related to the input signal vector by theequation:$y_{k} = {x_{k} - {\frac{1}{M}{\sum\limits_{j = 0}^{M - 1}x_{k - {j\; \frac{N}{M}}}}}}$where y_(k) represents an output sample of the FIR filter at index k,x_(k) represents a sample of the input signal at index k, N represents asample size of the input signal vector, and M is an integer.
 14. Themethod of claim 14, wherein M is selected to notch an interfering signalhaving a fundamental frequency and harmonic components from the inputsignal.
 15. The method of claim 13, wherein the linear FIR filter thatemulates a product of an inverse Fast Fourier Transform (IFFT) matrixand a selectively zeroed Fast Fourier Transform (FFT) matrix.
 16. Themethod of claim 13, wherein the selectively zeroed FFT matrix has selectrows set to zero.
 17. The method of claim 16, wherein the select rowsare selected based on an interfering signal.
 18. The method of claim 13,wherein the selectively zeroed FFT matrix has every M^(th) row set tozero.
 19. The method of claim 13, wherein N/M is an integer.